{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "UrWJ4QBTJRsv"
   },
   "source": [
    "# Working with Celestial Coordinates in WCS 1: Specifying, reading, and plotting\n",
    "\n",
    "## Authors\n",
    "\n",
    "Kris Stern, Kelle Cruz, Lia Corrales, David Shupe, Adrian Price-Whelan\n",
    "\n",
    "## Learning Goals\n",
    "\n",
    "1. Demonstrate two ways to build a `astropy.wcs.WCS` object\n",
    "2. Show an image of the Helix nebula with RA and Dec labeled\n",
    "3. Plot a scale bar on an image with WCS information\n",
    "\n",
    "## Keywords\n",
    "\n",
    "WCS, coordinates, matplotlib\n",
    "\n",
    "## Companion Contents\n",
    "1. \"An Introduction to Modern Astrophysics\" ([Carroll & Ostlie](https://ui.adsabs.harvard.edu/abs/2006ima..book.....C/abstract))\n",
    "2. [FITS WCS page at GSFC](https://fits.gsfc.nasa.gov/fits_wcs.html)\n",
    "\n",
    "## Summary\n",
    "\n",
    "This tutorial series aims to show how the content of Chapter 1 of \"An Introduction to Modern Astrophysics\" by Carroll and Ostlie can be applied to real life astrophysics research situations, using tools in the Astropy ecosystem. We will introduce two different approaches to building a `astropy.wcs.WCS` object, which contains meta-data that (in this case) defines a mapping between image coordinates and sky coordinates. The `astropy.wcs` subpackage conforms to the standards of the FITS World Coordinate System (WCS) used extensively by the astronomy research community. We will created a 2D WCS for an image of the iconic the Helix nebula (a planetary nebula) and display an image of the nebula with sky coordinates (here, equatorial, ICRS RA and Dec.) labeled. Finally, we will over-plot a scale bar on the Helix nebula image using WCS to give the reader a sense of the angular size of the image."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "ApMUBsVVJRsw"
   },
   "outputs": [],
   "source": [
    "from astropy.wcs import WCS\n",
    "from astropy.io import fits\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "AbDZ0uWTJRsz"
   },
   "source": [
    "## Section 1: Two ways to create an `astropy.wcs.WCS` object\n",
    "\n",
    "*World coordinates* serve to locate a measurement in some multi-dimensional parameter space. A World Coordinate System (WCS) specifies the physical, or world, coordinates to be attached to each pixel or voxel of an N-dimensional image or array. An [elaborate set of standards and conventions](https://fits.gsfc.nasa.gov/fits_wcs.html) have been developed for the Flexible Image Transport System (FITS) format ([Wells et al. 1981](https://ui.adsabs.harvard.edu/abs/1981A&AS...44..363W/abstract)). A typical WCS example is to specify the Right Ascension (RA) and Declination (Dec) on the sky associated with a given the pixel or spaxel location in a 2-dimensional celestial image ([Greisen & Calabretta 2002](https://ui.adsabs.harvard.edu/abs/2002A&A...395.1061G/abstract); [Calabretta and Greisen 2002](https://ui.adsabs.harvard.edu/abs/2002A&A...395.1077C/abstract)).\n",
    "\n",
    "The [`astropy.wcs` subpackage](https://docs.astropy.org/en/stable/wcs/) implements FITS standards and conventions for World Coordinate Systems. Using the `astropy.wcs.WCS` object and `matplotlib`, we can generate images of the sky that have axes labeled with coordinates such as right ascension (RA) and declination (Dec). This requires selecting the proper projections for `matplotlib` and providing an `astropy.visualization.WCSAxes` object.\n",
    "\n",
    "There are two main ways to initialize a `WCS` object: with a Python dictionary (or dictionary-like object, like a FITS file header) or with Python lists. In this set of examples, we will initialize an `astropy.wcs.WCS` object with two dimensions, as would be needed to represent an image.\n",
    "\n",
    "The WCS standard defines a set of keywords that are used to represent the world coordinate system for a given set of data (e.g., image). Here is a list of the essential WCS keywords and their uses; In each case, the integer $n$ denotes the dimensional axis (starting with 1) to which the keyword is being applied. In our examples below, we will have two image dimensions (axes), so $n$ will either be 1 or 2.\n",
    "\n",
    "* **CRVALn**: the coordinate value at a reference point (e.g., RA and DEC value in degrees)\n",
    "* **CRPIXn**: the pixel location of the reference point (e.g., CRPIX1=1, CRPIX2=1 describes the center of a corner pixel)\n",
    "* **CDELTn**: the coordinate increment at the reference point (e.g., the difference in 'RA' value from the reference pixel to its neighbor along the RA axis)\n",
    "* **CTYPEn**: an 8-character string describing the axis type (e.g., 'RA---TAN' and 'DEC---TAN' describe the typical tangent-plane sky projection that astronomers use)\n",
    "* **CUNITn**: a string describing the unit for each axis (if not specified, the default unit is degrees.)\n",
    "* **NAXISn**: an integer defining the number of pixels in each axis\n",
    "\n",
    "Some good references of the WCS standard can be found [here](https://fits.gsfc.nasa.gov/fits_wcs.html)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "A_uD9fzbUlet"
   },
   "source": [
    "### Method 1: Building a WCS object with a dictionary\n",
    "\n",
    "One way to define an Astropy `WCS` object is to construct a dictionary containing all of the essential information (i.e., specifying values for the WCS keywords listed above) that map the pixel coordinate space to the world coordinate space. \n",
    "\n",
    "In this example, we define two coordinate axes with:\n",
    "* A Gnomonic (tangent-plane) projection, which corresponds to the RA/Dec coordinate system\n",
    "* A reference location of (RA,DEC) = (337.52, -20.83), as defined by the **CRVALn** keys\n",
    "* The pixel at coordinate value (1,1) as the reference location (**CRPIXn** keys)\n",
    "* Units of degrees (**CUNITn = 'deg'**)\n",
    "* Pixel sizes of 1 x 1 arcsec (**CDELTn = 0.002778** in degrees)\n",
    "* An image size of 1024 x 1024 pixels (**NAXISn** key)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "3fRBbYmbJRs0"
   },
   "outputs": [],
   "source": [
    "wcs_input_dict = {\n",
    "    'CTYPE1': 'RA---TAN', \n",
    "    'CUNIT1': 'deg', \n",
    "    'CDELT1': -0.0002777777778, \n",
    "    'CRPIX1': 1, \n",
    "    'CRVAL1': 337.5202808, \n",
    "    'NAXIS1': 1024,\n",
    "    'CTYPE2': 'DEC--TAN', \n",
    "    'CUNIT2': 'deg', \n",
    "    'CDELT2': 0.0002777777778, \n",
    "    'CRPIX2': 1, \n",
    "    'CRVAL2': -20.833333059999998, \n",
    "    'NAXIS2': 1024\n",
    "}\n",
    "wcs_helix_dict = WCS(wcs_input_dict)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "AA4wMjYqJRs2"
   },
   "source": [
    "Now let's print the `WCS` object defined with a Python dictionary to check its content: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 194
    },
    "colab_type": "code",
    "id": "_3_wPImmJRs3",
    "outputId": "7629004b-9897-4c94-e681-1217fa8dcfb5"
   },
   "outputs": [],
   "source": [
    "wcs_helix_dict # To check output"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this demonstration (and below), we assume that we know all of the relevant WCS keyword values to specify. Typically, however, we will rely on software to produce these values for us. For example, WCS information is most often included automatically in FITS files produced by software used to take images with most instruments at astronomical observatories. In cases when the WCS information is provided for us in a FITS file, it is typically included in a FITS header, which, when read into Python, acts like a dictionary object. We demonstrate this later on in this tutorial."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "YlwK4xHQU6Un"
   },
   "source": [
    "### Method 2: Create an empty WCS object before assigning values\n",
    "\n",
    "Alternatively, we could initialize the `astropy.wcs.WCS` object, and assign the keyword values with lists corresponding to each respective axis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "bGI2DpWdJRs6"
   },
   "outputs": [],
   "source": [
    "wcs_helix_list = WCS(naxis=2)\n",
    "wcs_helix_list.wcs.crpix = [1, 1]\n",
    "wcs_helix_list.wcs.crval = [337.5202808, -20.833333059999998]\n",
    "wcs_helix_list.wcs.cunit = [\"deg\", \"deg\"]\n",
    "wcs_helix_list.wcs.ctype = [\"RA---TAN\", \"DEC--TAN\"]\n",
    "wcs_helix_list.wcs.cdelt = [-0.0002777777778, 0.0002777777778]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "ZjCbtC-bAtMr"
   },
   "source": [
    "Let's print the `WCS` object one more time to check how our values were assigned."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 194
    },
    "colab_type": "code",
    "id": "iwx8iECxJRs8",
    "outputId": "d8a4818f-13b0-433f-d17d-f660e207ba4c"
   },
   "outputs": [],
   "source": [
    "wcs_helix_list # To check output"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "rL7-Qut6AGDM"
   },
   "source": [
    "Note that when we initialize the WCS object this way, the `NAXIS` values are set to 0. To assign coordinates to our image, we will need to fix the shape of the `WCS` object array so that it matches our image. We can do this by assigning a value to the `array_shape` attribute of the `WCS` object:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "6PnSCW0_-pQs"
   },
   "outputs": [],
   "source": [
    "wcs_helix_list.array_shape = [1024, 1024]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "AszZGXxNBDoz"
   },
   "source": [
    "Now when we print the `WCS` object, we can see that the `NAXIS` values have been updated from the default size of 0 to 1024."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 194
    },
    "colab_type": "code",
    "id": "qdg6jLCa_Sck",
    "outputId": "8285ba5c-10e8-457f-bded-15c3faff782c"
   },
   "outputs": [],
   "source": [
    "wcs_helix_list"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "qKvxquTdEOo5"
   },
   "source": [
    "## Section 2: Show an image of the Helix nebula with RA and Dec labeled\n",
    "\n",
    "Most of the time we can obtain the required `astropy.wcs.WCS` object from the header of the FITS file from a telescope or astronomical database. This process is described below.\n",
    "\n",
    "### Step 1: Read in the FITS file\n",
    "\n",
    "We will read the FITS file containing an image of the Helix nebula from the `astropy-data` GitHub repository using the `astropy.io.fits` subpackage. The `astropy.io.fits.open()` function will load the contents of a FITS file into Python, and it accepts either a local file path or a URL (as is demonstrated here). This image (FITS file) was originally accessed from the [Digitized Sky Survey](https://archive.eso.org/dss/dss) but is provided in the `astropy-data` repository for convenience:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "MRVbLwvcE__Y"
   },
   "outputs": [],
   "source": [
    "header_data_unit_list = fits.open('https://github.com/astropy/astropy-data/raw/6d92878d18e970ce6497b70a9253f65c925978bf/tutorials/celestial-coords1/tailored_dss.22.29.38.50-20.50.13_60arcmin.fits')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "tTKAPz2pJzH0"
   },
   "source": [
    "FITS files are a binary file format that is mainly used by astronomers and can contain information arranged in many \"extensions,\" which contain both header information (e.g., metadata) and data (e.g., image data). We can check how many extensions there are in a FITS file, as well as view a summary of the contents in each extension, by printing the FITS object information."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 70
    },
    "colab_type": "code",
    "id": "n5T4-QCzJVKa",
    "outputId": "f7cb2735-6c58-44b3-a965-bd297e2bfc10"
   },
   "outputs": [],
   "source": [
    "header_data_unit_list.info()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "we-RjJ8PJo1w"
   },
   "source": [
    "This shows us that our FITS file contains only one extension, labeled 'PRIMARY' (or extension number 0). We will copy the image data from this extension to the variable `image`, and the header data to the variable `header`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "w4W_HiR-Ji6g"
   },
   "outputs": [],
   "source": [
    "image = header_data_unit_list[0].data\n",
    "header = header_data_unit_list[0].header"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "EUVTT2WdKCGd"
   },
   "source": [
    "We can print the FITS image header to screen so that all of its contents can be checked or utilized. Note that the WCS information for this information can be found near the bottom of the printed header, below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 1000
    },
    "colab_type": "code",
    "id": "XcerK_cZJ30W",
    "outputId": "f609c0ae-0116-45c6-fd89-921a208c9434"
   },
   "outputs": [],
   "source": [
    "header"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "g3FyYZ-bQ5b7"
   },
   "source": [
    "Please note that the *original* header (as downloaded from the DSS) violates the FITS WCS standards (because it includes both CDELTn keywords and a matrix of CD values; including deprecated PC-matrix keywords). The header has been cleaned up to conform to the existing standards."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "FzAtb_LWLyWn"
   },
   "source": [
    "### Step 2: Read in the FITS image coordinate system with astropy.wcs.WCS\n",
    "\n",
    "Because the header contains WCS information and acts like a Python dictionary, an Astropy `WCS` object can be created directly from the FITS header."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "3Wl5J0HBLoIp"
   },
   "outputs": [],
   "source": [
    "wcs_helix = WCS(header)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "Akv472vTZFLp"
   },
   "source": [
    "Let's print the `WCS` object to see what values were drawn from the header."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 194
    },
    "colab_type": "code",
    "id": "1C7y709IZNvl",
    "outputId": "6aee1d99-9f6f-4990-ef55-aa601c33ff0d"
   },
   "outputs": [],
   "source": [
    "wcs_helix"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "O10H9rmtQEnM"
   },
   "source": [
    "### Step 3: Plot the Helix nebula with sky coordinate axes (RA and Dec)\n",
    "\n",
    "The image data, `image`, is a 2D array of values, and by itself contains no information about the sky coordinates of the pixels. So, if we plotted the image by itself, the plot axes would show pixel values. (We will be using the `matplotlib` library for the plotting.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(10, 10))\n",
    "plt.imshow(image, origin='lower', cmap='cividis')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "All of the information that maps from these pixel values to sky coordinates comes from the WCS metadata, which we loaded into the `wcs_helix` object (from the FITS file header). This `WCS` object is built so that it can be provided to `matplotlib` with the  `projection` keyword, as shown in the call to `matplotlib.pyplot.subplot` below, in order to produce axes that show sky coordinate information instead of pixel values. We will also overlay a coordinate grid in ICRS equatorial coordinates by passing the sky coordinate frame name (here, \"icrs\") to the `ax.get_coords_overlay()` method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 607
    },
    "colab_type": "code",
    "id": "fnv48KDbOkxB",
    "outputId": "05775d9a-6663-4b01-a3d3-0d7c9440ea64"
   },
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(10, 10))\n",
    "ax = plt.subplot(projection=wcs_helix)\n",
    "plt.imshow(image, origin='lower', cmap='cividis', aspect='equal')\n",
    "plt.xlabel(r'RA')\n",
    "plt.ylabel(r'Dec')\n",
    "\n",
    "overlay = ax.get_coords_overlay('icrs')\n",
    "overlay.grid(color='white', ls='dotted')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "DVVT552iQYn4"
   },
   "source": [
    "## Exercise\n",
    "\n",
    "Copy the code block above and instead overlay a coordinate grid in Galactic coordinates."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 610
    },
    "colab_type": "code",
    "id": "XPA-2LOLQLnB",
    "outputId": "7eb7e024-b8a7-4d16-b0c3-0a1aa7c71d39"
   },
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(10, 10))\n",
    "ax = plt.subplot(projection=wcs_helix)\n",
    "plt.imshow(image, origin='lower', cmap='cividis', aspect='equal')\n",
    "plt.xlabel(r'RA')\n",
    "plt.ylabel(r'Dec')\n",
    "\n",
    "overlay = ax.get_coords_overlay('galactic')\n",
    "overlay.grid(color='white', ls='dotted')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "FexGrb2hXj85"
   },
   "source": [
    "## Section 3: Plot a scale marker on an image with WCS\n",
    "\n",
    "To add a scale marker (i.e., a line of a particular *angular* size) to the image of the Helix nebula, we will use the matplotlib `Axes.arrow` method to draw a line. \n",
    "\n",
    "First, we need to decide where to place the scale bar. In the example below, we define the center of the scale marker to be at `(RA, Dec) = (337 deg, -21.2 deg)`. \n",
    "\n",
    "We then use the `transform` attribute of `Axes.arrow` to draw our scale bars in degrees (instead of pixel coordinates). In this case, we draw a scale marker with a length of 0.1 degrees. The arrow method inputs are `ax.arrow(x, y, dx, dy, **kwargs)`, with `x` and `y` being the `RA` and `Dec` of the beginning of the line. We use `dx=0` so that there is no horizontal component in the bar, and `dy=0.1`, which gives the length of the arrow in the vertical direction. To ensure that the arrow is drawn in the J2000 ICRS coordinate frame, we pass `ax.get_transform('icrs')` to the `transform` keyword.\n",
    "\n",
    "Finally, we use `matplotlib.pyplot.text` to mark the length of the scale marker."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 607
    },
    "colab_type": "code",
    "id": "OkH12BWESpGp",
    "outputId": "adcffc9d-85fb-45f2-e38a-11f18db320d7"
   },
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(10, 10), frameon=False)\n",
    "ax = plt.subplot(projection=wcs_helix)\n",
    "ax.arrow(337, -21.2, 0, 0.1, \n",
    "         head_width=0, head_length=0, \n",
    "         fc='white', ec='white', width=0.003, \n",
    "         transform=ax.get_transform('icrs'))\n",
    "plt.text(337.05, -21.18, '0.1 deg', \n",
    "         color='white', rotation=90, \n",
    "         transform=ax.get_transform('icrs'))\n",
    "plt.imshow(image, origin='lower', cmap='cividis', aspect='equal')\n",
    "plt.xlabel(r'RA')\n",
    "plt.ylabel(r'Dec')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "N_CWLccBPXvu"
   },
   "source": [
    "## Exercise\n",
    "\n",
    "Make a horizontal bar with the same length. Keep in mind that 1 hour angle = 15 degrees."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "bIePbBar8uFv"
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.3"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
